library(foreign)
achrcourtjud <- read.dta(BayesianIACtHRMasterData.dta)
summary(achrcourtjud)

#set working directory
setwd("")
library(R2jags)

## Now define the vectors of the data matrix for JAGS:

achrJUD <- achrcourtjud$numcasejud_1
achrviol <- achrcourtjud$iacthr_numcase_1
jud_var <- achrcourtjud$ljivar_1
jud <- achrcourtjud$lji_1
fdi <- achrcourtjud$lfdi_1
elections <- achrcourtjud$v2xel_frefair_1
veto <- achrcourtjud$polconiii_1
exec <- achrcourtjud$xrcomp_1
speech <- achrcourtjud$speech_1
cs <- achrcourtjud$cs_1
nhri <- achrcourtjud$nhri_1
regionalcount <- achrcourtjud$regionalcountbin_1
gdp <- achrcourtjud$gdpcaplog
popmillions <- achrcourtjud$popmillionslog
cwar <- achrcourtjud$cwar
physint <- achrcourtjud$repressionlatent
physint_var <- achrcourtjud$repressionlatentvar
ccode <- achrcourtjud$ccode
id <- achrcourtjud$id
year <- achrcourtjud$year

N <- length(achrcourtjud$ccode)

## Read in the Courts data for JAGS

achrcourtjud.data  <- list("id", "achrJUD", "achrviol", "jud", "jud_var", "fdi", "elections", "veto", "exec", "speech", "cs", "nhri", "regionalcount", "physint_var", "gdp", "popmillions", "cwar", "physint", "N")

## Name the JAGS parameters

achrcourtjud.params <- c("alpha", "b1", "b2", "b3", "b4", "b5", "b6", "b7", "b8", "b9", "b1.0", "b1.1", "b1.2", "b1.3", "b1.4", "tau")

#Specify intial values (taken from linear regression in Stata)

achrcourtjud.inits1 <- list("alpha[id]"= -1, "b1"= -1, "b2"= -1, "b3"= -1, "b4"= -1, "b5" = -1, "b6"=-1, "b7"= -1, "b8"= -1, "b9" = -1, "b1.0" = -1, "b1.1" = -1, "b1.2" = -1, "b1.3" = -1, "b1.4" = -1)
achrcourtjud.inits2 <- list("alpha[id]"= 1, "b1"= 1, "b2"= 1, "b3"= 1, "b4"= 1, "b5" = 1, "b6"=1, "b7"= 1, "b8"= 1, "b9" = 1, "b1.0" = 1, "b1.1" = 1, "b1.2" = 1, "b1.3" = 1, "b1.4" = 1)
achrcourtjud.inits <- list(achrcourtjud.inits1, achrcourtjud.inits2)

achrjudfit <- jags(data=achrcourtjud.data, inits=achrcourtjud.inits,
                                parameters=achrcourtjud.params, n.chains=2, n.iter=100000, n.burnin=20000,
                                model.file="BayesIACtHRjudiciary.jags")

achrjudfit.upd <- update(achrjudfit, n.iter=200000)



###Diagnostics####
library(mcmcplots)
achrjudfit.mcmc <- as.mcmc(achrjudfit.upd)

traplot(achrjudfit.mcmc, parms = c("b1", "b2", "b3"), style = "plain", main = NULL, col = c("red", "black"))
denplot(achrjudfit.mcmc, parms = c("b1", "b2", "b3"), style = "plain", main = NULL, col = c("red", "black"))

caterplot(achrjudfit.mcmc, parms = c("b1", "b2", "b3", "b4", "b5", "b6", "b7", "b8", "b9", "b1.0", "b1.1", "b1.2", "b1.3", "b1.4", "b1.5"), labels = c("IACtHR*Jud", "IACtHR", "Jud",  "Aid", "ExecElec", "Access", "Checks", "Speech", "CS", "NHRI", "Democracy", "GDP (logged)", "Population (logged)", "CWar", "Regional"),
          title("", xlab = "Parameter Estimates"), quantiles = list(outer=c(0.05, 0.95), inner=c(0.05, 0.95)), style = "plain", abline(v = 0, lty = 2), pch = 19)

jagachrjud <- print(achrjudfit.upd, intervals=c(0.025, 0.05, 0.95, 0.975))
